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Abstract 

The use of exactly-solvable Richardson-Gaudin (R-G) models to describe the 
physics of systems with strong pair correlations is reviewed. We begin with a brief 
discussion of Richardson's early work, which demonstrated the exact solvability of 
the pure pairing model, and then show how that work has evolved recently into 
a much richer class of exactly-solvable models. We then show how the Richard- 
son solution leads naturally to an exact analogy between such quantum models and 
classical electrostatic problems in two dimensions. This is then used to demonstrate 
formally how BCS theory emerges as the large- TV limit of the pure pairing Hamilto- 
nian and is followed by several applications to problems of relevance to condensed 
matter physics, nuclear physics and the physics of confined systems. Some of the 
interesting effects that are discussed in the context of these exactly-solvable mod- 
els include: (1) the crossover from superconductivity to a fluctuation-dominated 
regime in small metallic grains, (2) the role of the nucleon Pauli principle in sup- 
pressing the effects of high spin bosons in interacting boson models of nuclei, and 
(3) the possibility of fragmentation in confined boson systems. Interesting insight 
is also provided into the origin of the superconducting phase transition both in 
two-dimensional electronic systems and in atomic nuclei, based on the electrostatic 
image of the corresponding exactly-solvable quantum pairing models. 
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I. INTRODUCTION 

Exactly-solvable models have played a major role in helping to elucidate the physics 
of strongly correlated quantum systems. Examples of their extraordinary success can be 
found throughout the fields of condensed matter physics and nuclear physics. In condensed 
matter, the most important exactly- solvable models have been developed in the context of 
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one-dimensional (ID) systems, where they can be classified into three families (Ha, 
The first family began with Bethe's exact solution of the Heisenberg model (JBethe, 



1996). 



1931). 



Since then a wide variety of ID models have been solved usin g the Bethe ansatz. A seconc 



1963; 



Tomonaea . 



famil y consists of the so-called Tomonaga-Luttinger models f)Luttinger , 
195oL which are solved by bosonization techniques and which have played an important 



role in revealing the non- Fermi-liquid propertie s of ferm i on sy stems in ID. For this reason, 



by 


Caloex 


TO 


( 


1962 


) and 


Sutherland 


( Haldane . 


1988; 


Sha.strv. 


M988 


), are 



applied to a variety of important problems, including the physics of spin systems, the quan- 
tum Hall effect, random matrix theory, electrons in ID, etc. 

Exactly-solvable models have also been developed in the field of nuclear physics, but from 
a different perspective. In these models, the Hamiltonian is written as a linear combination 
of the Casimir operators of a group decomposition chain chosen to represent the physics 



of a particular nuclear phase. One example is the SU(3) model of 



Elliott 



M958I ). which 



describes nuclear deformation and the associated rotational motion. Others include the 



three dynamical symmetry limits of the U(6) Interacting Boson Model (jlachello and Arima . 



1980T ). which respectively describe rotational nuclei (the SU(3) limit), vibrational nuclei (the 



U(5)limit) and gamma- unstable nuclei (the 0(6) limit). These models, all exactly solvable, 
have been extremely useful in providing benchmarks for a description of the complicated 
collective phenomena that arise in nuclear systems. 

Superconductivity is a phenomenon that is common to both nuclear physics and con- 
densed matter systems. It is typically descri bed by assuming a p airing Hamiltonian and 



treating it at the level of BCS approximation ((Bardeen et al 



1957), an approximation that 



explicitly violates particle number conservation. While this limitation of the BCS approx- 
imation has a negligible effect for macroscopic systems, it can lead to significant errors 
when dealing with small or ultrasmall systems. Since the fluctuations of the particle num- 
ber in BCS are of the order of y/~N (N being the number of particles), improvements of 
the BCS theory are required for systems w i th N ~ 100 particles. The number-projected 



BCS (PBCS) approximation ((Dietrich et al 



19641 ) has been developed and used in nuclear 



physics for a lo ng time and more recently ha s been applied in studies of ultrasmall supercon- 



ducting grains ((Braun and von Delft . 



19981 ). In the latter, it has been shown necessary to 
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go be yond PBCS approximation and to resort to the exact solution ([Dukelskv and Sierra . 



19991 ) to properly describe the crossover from the superconducting regime to the pairing 

solutio n of the pairing 



fluctuation regime. It was in this context that the exa ct numerica 



mode l, pu blished in a series of pap e r in t he sixties by 



1968) and 



Richardson 



1965 



1966a. 



Richardson and Sherman! ()196J) and independently disc ussed by Gaudin ( 1995), 



was rediscove red and applied successfully to small metallic grains ([Sierra et all 120001 ) . [For 



a review, see 



von Delft and Ralph! ( 2001 ).] The exact Richardson solution, though present 



in the nuclear physics literature since the sixties, was scarcely used until very recently , with 



but a few exceptions ([Bang and Krumlimde , 



1970; 



Hasegawa and Tasaki. 



mm 



1993). 



Soon after the initial application of the exact solution of the pairing model to ultra- 
small superconducting grains, it became clear that there is an intimate connection between 
Richards on's solu t ion an d a different family of exactly-solvable models known as the Gaudin 



magnet ((GaudinL I1976T). The connection came throu g 



the pairing model by 



i the proof of the integrability of 



Gambiaggio. Rivas and Saracenol (jl997l ). who identified the complete 



set of commuting operators of the model, the quantum invariants whose eigenvalues are 
the constants of motion. This made it possible to recast the pairing Hamiltonian as a 
linear combination of the quantum invariants. Furthermore, by establishing a relation be- 
tween the constants of motion of the pairing model and those of the Gaudin magnet, it 
became possible to generalize the pairing model to three classes of pairing-like models, all 
of which were i ntegrable and all o f which could be solve d exactly for both fermion and 



boson systems (jAmico et al 



2001 



Dukelskv et al 



20011 ). The potential importance of 



these exactly-solvable mode 



s in hig h T r su perconductivity and other fields of physics has 



recently been pointed out ((Heritieri 120011 ). It has also recently been show n how to ob 



tain the exact so 



2001 



Links et al. 



utions of these models using the Algebraic Bethe ansatz (jAmico et al 



2003 



von Delft and Poghossian . 



2002: 



Zhou et al. 



20021 ) . and a c onnec 



tion with Conformal Fie 
and 



Asorev et al 



d Theory and Chern-Simons theory has been established by 

(2nd). 



Sierra 



Fo 



Richardson 



lowin g their i n itial d iscussion of these exactly-solvable quantum models 
(|l977f) and iGaudinl ([19951 ) proposed an exact mapping between these models and a two- 
dimensional classical electrostatic problem. By exploiting this analogy, they were able to 
derive the thermodynamic limit of the exact solution, demonstrating that it corresponds 
precisely to the BCS solution. Recently, the same electrostatic mapping has been used 
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to provide a pictorial view of the transition to superconductivity in finite nuclei and t o 



suggest an alternative geometrical characterization of this transition (D ukelskv et al 



2002L 



imit based on the e lectrostatic image has 



2002) 



Furthermore, the validity of the thermodynamic 
been numerically checked for very large systems ((Roman et al. 

In this Colloquium, we review the recent progress that has been made in the use of 
exactly-solvable pairing and pairing-like models for the description of strongly-correlated 
quantum many-body systems. We begin with a review of the early work of Richardson and 
Gaudin, who first showed the exact solvability of such models. We then discuss the extension 
to a wider class of exactly-solvable models, building on the ideas arising from their integra- 
bility. We then discuss several recent ap plications, to ultrasmall superconducting grains 



( Schechter et al 



2001 



Sierra et al. 



ture (Dukelskv and Pittell 



20001 ) . to interacting boson models of nuclear struc 



20011 ). to e 



20011 ) , and to confined Bose systems ([Dukelskv and Schuckl . 120011 ) . 



ectrons in a two-dimensional lattice ( Dukelskv et al. 



II. THE MODELS OF RICHARDSON AND GAUDIN AND THEIR GENERAL- 
IZATION 

A. Richardson's exact solution of the pairing model 



The pairing interaction is the part of the fermion Hamiltonian responsible for the super- 
conducting phase in metals and in nuclear matter or neutron stars. It is also responsible for 
the strong pairing correlations in the corresponding finite systems, namely ultrasmall super- 
conducting grains and atomic nuclei. Its microscopic origin, however, depends on the system 
under discussion. In condensed matter, it derives from the exchange of phonons between 
the conduction electrons. In nuclear physics, it comes about due to the short-range nature 
of the effective nucleon-nucleon interaction in a nuclear medium and includes contributions 
both from the singlet-S and triplet-P channels of the nucleon-nucleon interaction. The main 
feature of the pairing interaction in both cases is that it co r relate s pairs of particles in time- 



reversed states. For recent reviews, see 



Sigrist and Uedal (jl99ll ) in condensed matter and 



Dean and Hjorth- Jensen (2003) in nuclear physics. 

Quite recently, the first direct sign of BCS superconductivity has been observed in a 
trapped degenerate gas of 40 K fermionic atoms. The pairing interaction in these trapped 
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dilute systems comes from the s-wave of the atom-atom interaction ([Heiselberg . 2003) and 
is characterized by the scattering length a. This property of the pairing interaction can 
be experimentally controlled by imposing an external magnetic field on the system, thereby 
tuning the location of the associated molecular Feshbach resonance. In this way, it is possible 
even to change the sign of the scattering length, which is a key to producing the observed 
fermionic superconductivity. 

We begin our discussion of Richardson's solution of the pairing model by assuming a 
system of N fermions moving in a set of L single-particle states /, each having a total degen- 
eracy Qi, and with an additional internal quantum number m that labels the states within 
the / subspace. If the quantum number / represents angular momentum, the degeneracy of 
a single-particle level I is Qi = 21 + 1 and —l<m<l. In general, however, I simply labels 
different quantum numbers. We will assume throughout this discussion of fermion systems 
that the ft/ are even so that for each state there is another obtained by time- reversal. The 
operators on which the pairing Hamiltonian is based are 



ni 



where a\ m (a, m ) creates (annihilates) a particle in the state (Im) and the state (7m) is the 
corresponding time- reversed state. The number operator fii, the pair creation operator A\ 
and the pair annihilation operator A[ close the commutation algebra 



41 = 25 w A}, \A h 41 = 26 w (ft, - 2nt) . (2) 



The corresponding algebra is SU{2). 

The most general pairing Hamiltonian can be written in terms of the three operators in 
Eq. (0) as 

H = Y,ein l + Y J V w A\A v . (3) 

i w 

Often a simplified Hamiltonian is considered, in which the pairing strengths V\y are replaced 
by a single constant g, giving rise to the pairing model (PM) or BCS Hamiltonian, 

ffp = £ e ifii + f£44' . (4) 
i A w 
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When g is positive, the interaction is repulsive; when it is negative, the interaction is at- 
tractive. 

The approximation leading to the PM Hamiltonian must be supplemented by a cutoff 
restricting the number of I states in the single-particle space. In condensed-matter problems 
this cutoff is naturally provided by the Debye frequency of the phonons. In nuclear physics, 
the choice of cutoff depends on the specific nucleus and on the set of active or valence orbits 
in which the pairing correlations develop. The cutoff in turn renormalizes th e strength of 



the e ffective pairing interaction that should be used within that active space ((Baldo et al 



1990). 



A generic state of M correlated fermion pairs and v unpaired particles can be written as 
|ni,n 2 r-',^,i/) = -=(4) ni (4) na ---(4) ni k) > ( 5 ) 



IN 

where M is a normalization constant. The number of pairs n\ in level I is constrained by the 
Pauli principle to be < 2n; + vi < O;, where V\ denotes the number of unpaired particles 
in that level. The unpaired state \v) = \ui, • • • Vl), with v = J2i Vu is defined such that 

Ai \v) = , ni \u) = vi \v) . (6) 

A state with v unpaired particles is said to have seniority v. The total number of collective 
(or Cooper) pairs is M = Y^i n i and the total number of particles is iV = 2M + v. 

The dimension of the Hamiltonian matrix in the Hilbert space of Eq. Q quickly exceeds 
the limits of large-scale diagonalization for a modest number of levels L and particles N. As 
an example, consider a problem involving L doubly-degenerate levels and M = N/2 pairs. 
One can readily carry out full diagonalization for up to L = 16 and M = 8, for which 
the full dimension of the Hamiltonian matrix is 12,870. For systems with up to L = 32 
and M = 16, the Lanczos method can still be used to obtain the lowest eigenvalues. But 
for larger values of L and M, such methods can no longer be used and it is necessary to 
find alternative methods to obtain physical solutions. BCS is an example of an alternative 
approximate method. Here we focus instead on finding the exact solution, but without 
numerical diagonalization. 

In spite of the apparent complexity of the problem, Richardson showed that the exact 
unnormalized eigenstates of the hamiltonian of Eq. Q can be written as 



S 



\^) = B\bI--B\ [ \u) , (7) 

where the collective pair operators B a have the form appropriate to the solution of the 
one-pair problem, 

In the one-pair problem, the quantities E a that enter Eq. (jHJ) are the eigenvalues of the PM 
Hamiltonian i.e., the pair energies. Richardson proposed to use the M pair energies E a 
in the many-body wave function of Eq. ((7|) as parameters which are then chosen to fulfill 
the eigenvalue equation H P = E 

We will not repeat here th e derivation of the set of equation s that the pair energies must 



fulfill, referring the reader to iRichardson and Sherman! ()1964f ). In Sect. IV, we will return 



to the eigenvalue problem for more general Hamiltonians, of which the pairing Hamiltonian 
is a particular case. 

The key conclusions from Richardson's derivation are as follows: 

• The state given in Eq. (JJJ) is an eigenstate of the PM Hamiltonian (jlj if the M pair 
energies E a satisfy the set of M nonlinear coupled equations (called the Richardson 
equations) 

i-4^E^V + ^ E ^V = o ; 0) 

, lei - z a moL) h a - tip 
where di = — is related to the effective pair degeneracy of single-particle level I. 

• The energy eigenvalue associated with a given solution for the pair energies is 



£ = E^ + E^- ( 10 ) 

l a 

Because the Richardson method reduces the problem to solving non-linear coupled equa- 
tions, it can be used for systems well beyond the limits of either exact numerical diagonal- 
ization or the Lanczos algorithm. For example, the method can be used to obtain exact 
solutions for systems with L = 1000 and M = 500, for which the dimension is 2.7 x 10 299 . 
While we cannot obtain all the solutions for such a problem, we can relatively easily obtain 
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the lowest few for any value of the coupling constant. This means that we can study the 
quantum phase transition for pairing using the Richardson algorithm but not any phase 
transitions associated with temperature. To do the latter, we would need to develop the 
Thermodynamic Bethe ansatz for application to the pairing model, as will be discussed 
further in the summary given in Sect. VII. 



B. The Gaudin magnet 



Inspired by Richardson's solution of the PM, and building on his previous work on the 



Bethe method for solving one-dimensional problems, iGaudinl (|1976( ) proposed a family of 
fully integrable and exactly-solvable spin models. The Gaudin models are based on the 
SU(2) algebra of the spin operators, 



K a , K p ] = 2ie a ^K~< , (11) 

where K a = a a (a = 1, 2, 3) are the Pauli matrices. 

A quantum model with L degrees of freedom is integrable if there exist L independent, 
global Hermitian operators that commute with one another. This condition guarantees the 
existence of a common basis of eigenstates for the L operators, called the quantum invariants, 
and for their eigenvalues, the constants of motion. 

Since the SU(2) algebra has one degree of freedom, the most general set of L Hermitian 
operators, quadratic in the spin operators and global in L spins, are 

Ri= E X>3 • (12) 

where the are 3L (L — 1) real coefficients. To define an integrable model, i.e., to satisfy 
the commutation relations [Hi, Hj] = 0, these coefficients must satisfy the system of algebraic 
equations 

K W ]k + w ji w ik ~ w ?kW% = • ( 13 ) 

Gaudin proposed two conditions to solve this system of equations. The first was anti- 
symmetry of the w coefficients, 
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K = • (14) 

The second was to express the w coefficients as an odd function of the difference between 
two real parameters, 

Kj = fa (Vi ~ Vj) , (15) 

in order to fulfill Eq. (JUJ). 

The most general solution of Eq. ()13|) subject to Eqs. (pPfj) and (|15|) can be written in 
terms of elliptic functions. We further restrict here to the case in which the total spin in 
the z direction, S z = \ J2i Kf, is conserved, i.e., the L operators Hi (Gaudin Hamiltonians) 
commute with S z . Conservation of S z requires that wfj = = and w-j = Yij, in terms 
of two new matrices X and Y. In such a case, the integrability conditions of Eq. (|T3|) reduce 
to 

YijXj k + Y ki Xj k + X ki Xij = . (16) 
There are three classes of solutions to Eq. (fTB|) : 



I. The rational model 



II. The trigonometric model 



Xij — Y^ — (17) 



X - = sm( Vl - Vj ) ' Y * = ( 18 ) 
III. The hyperbolic model 

X l3 = sinh( l_^ , Y tJ = coth ( Vi - Vj ) (19) 

We will hold off presenting details on the solution of the three families of Gaudin models 
until Sect. II. D, where we discuss the generalization of the Richardson- Gaudin (R-G) models. 
At that point we will see that the solutions of the Gaudin models are based on precisely the 
same ansatz as used by Richardson in his solution of the pairing model. 
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The most general integrable spin hamiltonian that can be written as a linear combination 
of the Gaudin integrals of motion in Eq. (JT2*Jl is 



E (C« - 0) fa [ K f K ! + K " K f\ + IW*; } • 



(20) 



The above Hamiltonian, which models a spin chain with long-range interactions, has a 
total of 2L free parameters. There are L rj's, which define the X's, and L ('s. It can 
be readily confirmed that the rational family gives rise to an XXX spin model while the 
trigonometric and hyperbolic families correspond to an XX Z spin model. To the best of our 
knowledge, there have been no physical applications of the Gaudin magnet, though there 
are indications that when the 2L free parameter s of the model are chosen at ra ndom the 
Gaudin magnet behaves as a quantum spin glass ((Arrachea and Rozenberd . 120011 ) . 

The Richardson model has a natural link to superconductivity and the Gaudin model to 
quantum magnetism, both very important concepts in contemporary physics. Despite these 
facts, neither model received much attention from the nuclear or condensed-matter commu- 
nities for many years. On the other hand, the Gaudin models h ave played an i mpor tant role 
in aspects of quantum integrability. For a recent reference, see lGould et al. ( 2002 ). 



C. Integrability of the pairing model 



Despite t he fact that the Richardson and Gaudin models are so similar, it was not until 



the work of 



Cambiaggio. Rivas and Saracend ([19971 ) (CRS) that a precise connection was 



established. We now discuss their work and show how it provides the necessary missing link. 

The key point of CRS was to show that the PM is integrable by finding the set of 
commuting Hermitian operators in terms of which the PM Hamiltonian could be expressed 
as a linear combination. Finding the complete set of common eigenvectors of those quantum 
invariant operators is then equivalent to finding the eigenvectors of the PM Hamiltonian. 

In their derivation, they began by introducing a pseudo-spin representation of the pair 
algebra, advancing a con nection between p airing phenomena and spin physics that had been 

( 19581) . 



pointed out long ago by 



Anderson 
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The elementary operators of the pair algebra, defined in terms of the generators of the 
SU (2) pseudo-spin algebra, are 



1 1 

K i = o a L a im ~ 7^ 



K7 



m Im 



(*r)' 



(21) 



The operator i^ 4 " creates a pair of fermions in time-reversed states. The degeneracy Qi 
of level I is related to a pseudo-spin 5/ for that level according to Qi = 2Si + 1. The three 
operators in Eq. (j2Tj) close the SU (2) commutation algebra, 



Kf, Kf \ = ±5 W K± , \K+, Kf] = 2S w Kf . 



(22) 



The ^C/ (2) group for a level I has one degree of freedom and its Casimir operator is 



(23) 



2 V / | 

For a problem involving L single-particle levels, there are obviously L degrees of freedom. 

Guided by previous work on the two-level PM, CRS considered the following set of oper- 
ators: 



'0 7^0 



i (^+^7 + KfK+ 



+ K°K$\ . (24) 

They showed that these operators are (1) Hermitian, (2) global, in the sense that they are 
independent of the Hilbert space, (3) independent, in the sense that no one can be expressed 
as a function of the others, and (4) commute with one another. Furthermore, there are 
obviously as many Ri operators as degrees of freedom. The set of L such operators thus 
fulfills the conditions required for the quantum invariants of a fully integrable model. Finally, 
they showed that the PM Hamiltonian of Eq. (@J) can be written as a linear combination of 
the Ri according to 



H P = 2Y,eiRi + C , (25) 
i 
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where C is an uninteresting constant. 

We now return to the relationship between the Richardson PM and the Gaudin models. 
This can be done by focussing on Gaudin's rational model and comparing its quantum 
invariants to those of CRS [see Eq. (J24)) ] for the pairing model. As can be readily seen, 
the two are very similar except that the quantum invariants of Gaudin's rational model 
are missing a one-body term or equivalently a linear term in the spin operators. As shown 
by CRS, this term preserves the commutability of the R operators and therefore generalizes 
Gaudin's rational model. In the following subsection, we discuss the solution of the so-called 
generalized Richardson- Gaudin models that emerge when this term is added. 

D. Generalized Richardson- Gaudin models 

The generalization of the Richardson and Gaudin models we now discuss proceeds along 
two distinct lines. First, we no longer limit our discussion to Gaudin's rational model, but 
now generalize to all three types (rational, trigonometric and hyberbolic). We will see that 
all three can be generalized by the inclusion of a linear term in the quantum invariants. 
Second, we generalize our discussion to include boson models as well as fermion models. 

We begin by considering the generalization to boson models. For boson systems, the 
elementary operators of the pair algebra are 

1 1 

^ m ^ 

K = ^E«L«L=(^r) t • (26) 

Note the similarity with the fermion pair operators of Eq. (j21|) . The only differences are that 
(1) there is a relative minus sign between the two terms in the operator K°, and (2) there 
is no longer a restriction to Qi being even. The latter point follows from the fact that for 
bosons a single-particle state can be its own time-reversal partner. As an example, consider 
scalar bosons confined to a 3D harmonic oscillator potential, as will be discussed further in 
Sect. VI. D. The degeneracy associated with a shell having principal quantum number n is 
&>n = {n+ 1) ■ (n+2)/2. The lowest two shells (n = and 1) have odd degeneracies, whereas 
the next two (n = 2 and 3) have even degeneracies 

The set of operators in Eq. (|2l)jl satisfy the commutation algebra 
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= -26 U 'Kf , (27) 

appropriate to SU(1,1). 

In subsequent considerations, we will often treat fermion and boson systems at the same 
time. To facilitate this, we can combine the relevant SU(2) commutation relations for 
fermions and SU(1,1) relations for bosons into the compact form 



± 



= &S U >K? . (28) 

When both types of systems are being treated together, we follow a convention whereby the 
upper sign refers to bosons and the lower sign to fermions. 

Following earlier discussion, we now consider the most general Hermitian and number- 
conserving operator with linear and quadratic terms, 



2 

TY U ,K°K$\ . (29) 



Note that this is a natural generalization of Eq. ()24j) . but now appropriate to both boson 
and fermion systems. 

We next look for the conditions that the matrices X and Y in Eq. ()29j) must satisfy 
in order that the R operators commute with one other. Surprisingly, the conditions are 
precisely those derived by Gaudin and presented in Eq. (jTBjl. despite the fact that the 
quantum invariants now include a linear term and that they now represent both boson and 
fermion systems. 

Here too three families of solutions derive, which can be written in compact form as 

sin [7 [rji -r)j)\ 

where 7 = corresponds to the rational model, 7 = 1 to the trigonometric model and 7 = 2 
to the hyperbolic model. The three limits are completely equivalent to those presented for 
the Gaudin magnet in Eqs. (jl7H19|) . 
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The next step is to find the exact eigenstates common to all L quantum invariants given 
inEq. flU, 

i2j|*) = ri|*) . (31) 

This can be accomplished by using an ansatz similar to the one used by Richardson [see Eq. 
(JZj)] to solve the PM, namely 

M L 

m=J[Bl\u) , Bl = Y,Ui{E a )Kt . (32) 

a=l i=l 

The wave function in Eq. (}3*2*j) is a product of collective pair operators B^, which are 
themselves linear combinations of the raising operators Kf that create pairs of particles in 
the various single-particle states. Note that they are analogous to the Richardson collective 
pair operators of Eq. (jSJ) • 

We now present explicitly the solutions to the eigenvalue equations for the three models. 
As we will see, it is possible to present the hyperbolic and trigonometric results in a single 
set of compact formulae, by using the notation sn for sin or sink, cs for cos or cosh and 
ct for cot or coth. These are not to be confused with elliptic functions. The rational model 
solutions are of a somewhat different structure and thus are presented separately. For each 
family, we first give the amplitudes that define the collective pairs in terms of the free 
parameters rji and the unknown pair energies E a , then the set of generalized Richardson 
equations that define the pair energies E a , and finally the eigenvalues of the quantum 
invariants B4. 

I. The rational model 



'^ErVrtE^-^ (34) 



n = di 



Vi Vj a ^Vi E a 



(35) 
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II and III. The trigonometric and hyperbolic models 



Ui (E a ) = — — — - , (36) 

sn {E a - 2r}i) 



lT4gJ2 d 3 ct{E a - 2tft) ±4g ct ( E P ~ E «) = °> ( 37 ) 



Ti = di {l =F 2g 



J2 dj ctfa-Vj) 



-2j2ct(E a -2r ]i ) 



(38) 



In Eqs. (j33H38j) . the quantity d\ is now given by 



* = i ± t ■ (39) 

Given a set of parameters r/, and a pairing strength g, the pair energies E a are obtained 
by solving a set of M coupled nonlinear equations, either those of Eq. (J3"4"j) for the rational 
model or those of Eq. (|37|) for the trigonometric or hyperbolic models. In the limit g — > 0, 
Eqs. (|34| IHTj) can only be satisfied for E a — > 2^j. In this limit, the corresponding pair 
amplitudes Ui(E a ) in Eqs. (JBlfl l3T)j) become diagonal and the states of Eq. (j3~2j) reduce to a 
product of uncorrelated pairs acting on an unpaired state. The ground state, for example, 
involves pairs filling the lowest possible states and no unpaired particles. Excitations involve 
either promoting pairs from the lowest paired states to higher ones or by breaking pairs and 
increasing the seniority. In this way, we can follow the trajectories of the pair energies E a 
that emerge from Eqs. ()34[ l3*7jl as a function of g for each state of the system. 

For boson systems the pair energies are always real, whereas for fermion systems the 
pair energies can either be real or can arise in complex conjugate pairs. In the latter case, 
there can arise singularities in the solution of Eqs. (}3"4"l IBTjl for some critical value of the 
pairing strength q r , when two or more pair energies acquire the same value. It was shown by 



Richardson! (jl965T ) that each of these critical g values is related to a single-particle level i and 
that at the critical point there are 1 — 2di pair energies degenerate at 2r)i. These singularities 
of course cancel in the calculations of energies, which do not show any discontinuity in the 
vicinity of the critical points. However, the numerical solution of the nonlinear equations 
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or (|37p may break down for values of g close to the singularities, making impractical 
the method of following the trajectories of the pair energ ies E n from the weak co upling limit 
to the desire value of the pairing strength g. Recently, (|Rombouts et all |2003) proposed a 
new method based on a change of variables that avoids the singularity problem opening the 
possibility of finding numerical solutions in the general case. 

The eigenvalues of the R operators, given by Eqs. (|3*5)l or P5]l. are always real since 
the pair energies are either real or come in complex conjugate pairs. Each solution of 
the nonlinear set of equations produces an eigenstate common to all Ri operators, and 
consequently to any Hamiltonian that is written as a linear combination of them. The 
corresponding Hamiltonian eigenvalue is the same linear combination of eigenvalues. 



III. THE ELECTROSTATIC MAPPING OF THE RICHARDSON-GAUDIN MOD- 
ELS 



We now introduce an exact mapping between the integrable R-G models and a classi- 
cal electrostatic p roble m in two dimen sions. Our derivation builds on the earlier work of 
Richardson! ()1977l ) and iGaudinl ()1995l ). who used this electrostatic analogy to show that 
the exact solution of the PM Hamiltonian agrees with the BCS approximation in the large- 
N limit. For simplicity, we concentrate here on the exact solution for the rational family 
(7 = 0). The exact solution of the trigonometric and the hyperbolic families can be reduced 
to a set of rational equations by a proper transform ation and then also interpreted as a 
classical two-dimensional problem (|Amico et all 12002). 

Let us assume that we have a two-dimensional (2D) classical system composed of a set 
M free point charges and another set of L fixed point charges. For reasons that will become 
clear when we use these electrostatic ideas as a means of studying quantum pairing problems, 
we will refer to the fixed point charges as orbitons and the free point charges as pairons. 

The Coulomb potential due to the presence of a unit charge at the origin is given by the 
solution of the Poisson equation 



V V (r) = -2tt5 (r) . (40) 

The Coulomb potential V (r) that emerges from this equation depends on the space 
dimensionality. For one, two and three dimensions, respectively, the solutions are 
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r ID 

V (r) ~ I In (r) 2D (41) 
1/r 3D 

In 3D, the Coulomb potential has the usual 1/r behavior, but in 2D it is logarithmic. In 
practice, there are no 2D electrostatic systems. However, there are cases of long parallel 
cylindrical conductors for which the end effects can be neglected so that the problem can be 
effectively reduced to a 2D plane. In our case, however, we will simply use the mathematical 
structure of the idealized 2D problem to obtain useful new insight into the physics of the 
quantum many-body pairing problem. 

For the purposes of this discussion, we map the 2D xy-plane into the complex plane by 
assigning to each point r a complex number z = x + iy. Let us now assume that the pairons 
have charges q a and positions z a , with a = 1, • • ■ , M, and that the orbitons have charges 
qi and positions Zi, with % = 1, • • • , L. Since they are confined to a 2D space, all charges 
interact with one another through a logarithmic potential. Let us also assume that there is 
an external uniform electric field present, with strength e and pointing along the real axis. 
The electrostatic energy of the system is therefore 



M L L M 

U = e Q a Re(z a ) + e ^ q j Re{z j ) ln 

a=l j=l j=l a=l 



- g 1<*QP 111 \ Z <* - Z P\ ' \ H Hlj 111 \ Z i ~ Z 3 I • ( 42 ) 

If we now look for the equilibrium position of the free pairons in the presence of the fixed 
orbitons, we obtain the extremum condition 



tE^-E^' («) 



It can be readily seen that the Richardson equations for the rational family (|34jl coincide 
with the equations for the equilibrium position of the pairons (|43J1 . if the. pairon charges are 
q a = 1, the pairon positions are z a = E a , the orbiton charges are equal to the effective 
level degeneracies di (which for the ground configuration are ±f^/4), the orbitons positions 
are Zi = 2r]i, and the electric field strength is e = ±^k 
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As a reminder, the T]iS are free parameters defining the integrals of motion from which 
the generalized R-G Hamiltonian is derived [see Eq. ([15)1] . In the case of the pure PM 
Hamiltonian, which is a specific Hamiltonian in the rational family, they are the single- 
particle energies of the active levels. In any subsequent discussion in which they are used, 
we will thus refer to them as effective single-particle energies. 

Table I summarizes the relationship between the quantum PM and the classical elec- 
trostatic problem implied by the above analogy. From the above discussion, we see that 
solving the Richardson equations for the pair energies E a is completely equivalent to finding 
the stationary solutions for the pairon positions in the analogous classical 2D electrostatic 
problem. 

TABLE I: Analogy between a quantum pairing problem and the corresponding 2D classical 

electrostatic problem. 



Quantum Pairing Model 


Classical 2D Electrostatic Picture 


Effective single particle energy rji 


Orbiton position Zi = 2r\i 


Effective orbital degeneracy di 


Orbiton charge qi = di 


Pair energy E a 


Pairon position z a = E a 


Pairing strength g 


Electric field strength e = ±j- 



Assuming that the real axis is vertical and the imaginary axis horizontal, and taking 
into account that the orbiton positions are given by the real single-particle energies, it is 
clear that they must lie on the vertical axis. The pairon positions are not of necessity 
constrained to the vertical axis, but rather must be reflection symmetric around it. This 
reflection symmetry property can be readily seen by performing complex conjugation on 
the electrostatic energy functional (|42|) . As a consequence, a pairon must either lie on the 
vertical axis (real pair energies) or must be part of a mirror pair (complex pair energies). 

The various stationary pairon configurations can be readily traced back to the weakly 
interacting system (g — > 0). In this limit, the pairons for a fermionic system are distributed 
around (and very near to) the orbitons, thereby forming compact artificial atoms. The num- 
ber of pairons surrounding orbiton i cannot exceed |2dj|, as allowed by the Pauli principle, 
i.e. we cannot accommodate in a single level more particles than its degeneracy permits. The 
lowest-energy (ground-state) configuration corresponds to distributing the pairons around 
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the lowest position orbitons consistent with the Pauli constraint. We then let the system 
evolve gradually with increasing g until we reach its physical value. Of course, the Pauli 
limitation does not apply to boson systems, which will also be discussed in the applications. 

To illustrate how the analogy applies for a specific quantum pairing problem involving 
fermions, we consider the atomic nucleus 114 Sn. This is a semi-magic nucleus, which can be 
modelled as 14 valence neutrons occupying the single-particle orbits of the N = 50 — 82 shell. 
Furthermore, it can be meaningfully treated in terms of a pure PM Hamiltonian with single- 
particle energies extracted from experiment. We will return to this problem briefly in Sect. 
VLB. For now we simply want to use this problem to illustrate the relationship between the 
quantum parameters and the electrostatic parameters for the analogous classical problem. 

In Table II, we list the relevant single-particle orbits of the 50 — 82 shell in the third 
column, their corresponding single-particle energies in the first column and the associated 
degeneracies in the second. Each level corresponds to an orbiton in the electrostatic problem, 
with the position of the orbiton given in the fourth column (at twice the single-particle energy 
of the corresponding single-particle level). Note that this is always pure real, meaning that 
in the 2D plot each orbiton has y = 0. Finally, in the fifth column we give the charge of the 
orbiton, which is simply related to the degeneracy of the level according to the prescription 
in Table I. 

TABLE II: Relationship between the quantum pairing problem for the nucleus UA Sn and 
the corresponding 2D classical problem. The notation "s.p." is shorthand for 
"single-particle". The s.p. energies are in MeV. 



s.p. energy 


s.p. degeneracy 


s.p. level/orb iton 


orbiton position 


orbiton charge 


0.0 


6 


d-5/2 


0.0 


-1.5 


0.22 


8 


97/2 


0.44 


-2.0 


1.90 


2 


Sl/2 


3.80 


-0.5 


2.20 


4 


d>3/2 


4.40 


-1.0 


2.80 


12 


hn/2 


5.60 


-3.0 



To illustrate the weak-pairing limit discussed above, we show in Fig. ^the pairon positions 
associated with the electrostatic solution for 114 Sn calculated for a pairing strength of g = 
—0.02 MeV, well below the strength at which the superconducting phase transition sets in. 
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Solid lines connect each pairon to its nearest neighbor. As we can readily see, the seven 
pairons in this case indeed distribute themselves very near to the lowest two orbitons, with 
three forming an artificial atom around the d§/2 and four forming an artificial atom around 
the 5f 7/2 . 



IV. THE LARGE-iV LIMIT 



We now discuss how the electrostatic mapping of the previous section can be used to study 
the exact solution of the pairing problem in the large- iV or thermodynamic limit. We focus 
on fer mion systems and consider the PM (or BCS) Hamiltonian used by 
f)1996l ) to describe the physics of ultrasmall superconducting grains, 



von Delft et al 



H B cs = 2 J2 e j a )a a jo ~ G Y1 

3 c=± 3 j' 



aj + aj_a,ji-aj'- 



(44) 



where aj± (a]±) are annihilation (creation) operators in the time-reversed single-particle 
states \j ±), both with energies Sj = €j/2, and G is the BCS coupling constant. Thus, ej 
denotes the energy of a pair occupying th e level j and €j ^ e,- for i ^ j. 



This model was solved analyticall y by 



Richardson and Sherman! (J1964J) and numerically 



up to L = 32 single-particle levels bv iRichardsonl (|1966af ). The seniority- zero eigenstates for 
a system of M fermions depend on a set of parameters E v {y — 1, . . . , M) (the pair energies) 
that are, in general, complex solutions of the M coupled algebraic Richardson equations 



1 

G 



M 



E,. 



- E 



E — E 



V = h 



M 



(45) 



The energy E associated with a given solution is given by the sum of the resulting pair 
energies [see Eq. (II OJ) ]. The ground state is given by the solution of Eq. (|43j) with the lowest 
value of E. 

In Fig. |2l we plot the solution of Eq. (J45|) for a model of equally-spaced energy levels, 
€j = d(2j — L — 1), j — 1, . . . , L, where d — uj/L is the single-particle level spacing and lj 
is twice the Debye energy. The calculation is done at half filling for M = L/2 = 8. [Note: 
At half filling, the number of levels L is equal to the number of particles N.] For small 
values of the coupling constant g = GL all the solutions E^ are real, but as we approach 
some critical value g Cj i the two real roots that are closest to the Fermi level coalesce at g c % 
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and for g > g Cj i develop into a complex conjugate pair. The same phenomenon happens to 
other roots as g is further increased, until ev entu ally all o f the r oots form complex pairs. 
This fact was observed by iRichardsonl ()1977l ) and iGaudinl ()1995l ) and suggested a way to 
analyze systems with a large number of particles, where the exact solution must converge 
asymptotically to the BCS solution. 

Figure |3] shows the solutions to Eq. (|45|) for a system with a much larger number of 
particles, M = N/2 = 100 pairs, and for three values of g. For g = 1.5, the roots form 
an arc which ends at the points 2A ± 2iA, where A is the BCS chemical potential and A 
the BCS gap. For g = 1.0 and 0.5, the set of roots consist of two pieces, one formed by an 
arc T with endpoints 2 A ± 2iA which touches the real axis at some point Ea, and a set of 
real roots along the segment [— uj,£a]- As g decreases, the latter segment gets progressively 
larger while the arc becomes smaller and eventually shrinks to a point when g = 0. 

The soli d lines in the figure are the results obtained from the algebraic equations derived 
by IGaudinl (jl995T ) in the large- N limit. Note that they are in excellent agreement with the 
results obtained by numerically solving Eq. (|45jl . 

We now show how to make the connection between the large- N limit of the PM problem 
and BCS theory more precise, by applying the electrostatic analogy introduced in Sect. IIII1 
We will consider the limit in which L — > oo, while keeping fixed the following quantities: 



r- 9 



M 

P = T 



(46) 



Assuming that the pair energies organize themselves into an arc T which is piecewise 
differentiate and symmetric under reflection on the real axis, Eq. ()45|) (the Richardson 
equation) in the continuum limit is converted into the integral equation 



POO de 



rjQ m 1 

2G 



0. 



(47) 



where p(e) is the energy density associated with the energy levels that lie in the interval 
Q = [—oo, u] and satisfies 



f L 

/ P( e ) rfe = 77 » 
Jn 2 

while r(£) is the density of the roots that lie in the arc V and satisfies 



(48) 
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jfr(£) \d£\ = M, 
J r £r(£M\ = E. 



(49) 
(50) 



The l ast eq uation is a consequence of Eq. (JTUJl . The solution of Eq. (}4*Tj) was given by 

Gaudin (1995) using techniques of complex analysis. We now summar ize his main results, 

I I I I 

which from a different perspective were also given bv iRichardsonl (|1977l ). A detailed deriva- 
tion of the continuum limit tog ether with a compari son with numerical results for large but 



finite systems was presented by 



Roman et al 



Introducing an "electric field", and studying its properties in the vicinity of the arc T, 
one can show that Eq. ()47j) yields the well-known BCS gap equation 



p(e) de 



A) 2 + A 2 



1 

G 



(51) 



that Eq. (f4*9~j) becomes the equation for the chemical potential 



M 



1 - 



- A 



A) 2 + A 2 



p(e)de, 



(52) 



and that Eq. (|50|) gives the BCS expression for the ground-state energy, 



E 



G + Jn[ 1 ^/(|-A) 2 + A 2 



p(e) e de. 



(53) 



Thus by using the electrostatic analogy for the quantum pairing problem we are able to 
demonstrate how the BCS equations emerge naturally in the large- N limit. 



V. THE ELEMENTARY EXCITATIONS OF THE BCS HAMILTONIAN 

Most of the studies to date of excited states of the BCS Hamiltonian [see Eq. (JHJ)] based 
on Richardson's exact solution have dealt with a subclass of these states, namely those 
obtained by breaking a single Cooper pair. In a case involving doubly-degenerate levels 
only, the levels in which the broken pair reside can no longer be occupied by the collective 
pairs. Those levels are thus blocked and this is reflected in the Richardson equations by them 
having effective degeneracy di = 0. For more general pairing problems, with degeneracies 
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larger than 2, a broken pair will not completely block a single-particle level, but rather will 
increase the seniority ui of the level and give rise to a reduced effective degeneracy di = y — 
There is also another class of collective excitations that arises without changing the seniority. 
These excitations, known in nuclear physics as pairing vibrations, correspond to the different 
solutions of the Richardson equations for a given seniority configuration z/j. 

In a mean-field treatment of the same Hamiltonian, the lowest such excitations of both 
types correspond to two quasi-particle states in BCS approximation, which are then mixed 
by the residual interaction in Random Phase Approximation (RPA). While this bears no 
obvious resemblance to the Richardson approach for these elementary excitations, it is clear 
that they should coincide in the large- N limit. 

To investigate this relation, we focus, for simplicity, on the BCS Hamiltonian of Eq. (|44j) 
and study systematically its excited states for systems with a fixed number of particles. Once 
we know the excited states, we can try to interpret them in terms of elementary excitations 
characterized by definite quantum numbers, statistics and dispersion relations. A generic 
excited state will them be given by a collection of elementary excitations. 

The elementary excitations within the exact Richardson approach are associated with the 
number of pair energies Ng, either real or in complex conjugate pairs, that in the large- q 



limit stay finite within t 



2003; 



r e inte rval between neighboring single-particle levels ([Roman et al 



Yuzbashvan et al 



2003). We display this behavior in Fig. 0] for the ground state 
(GS) and the first two excited states for a system with L = 40 single-particle levels at half 
filling (M = 20) as a function of g. Note that for the ground state all pair energies become 
complex and their real parts go to infinity as the coupling strength becomes infinite. Thus 
Ng = for this state. Correspondingly we find that for the first and second excited states 
the number of pair energies that remain finite are Nq = 1 and 2, respectively. 

In the large-g limit, the pair energies that stay finite (El) satisfy the Gaudin equation 



N G 

E-^J- E ^=0. v = l,...,No, (54) 



while the remaining M — Ng pair energies ( El) go to infinity and ar e the solution of the 



generalized Stieltjes problem encountered by Sh astrv and Dharl ([20011 ) in the study of the 



excitations of the ferromagnetic Heisenberg model, 
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I j m—Nq 2 

- + t^+ y tt- ff = 0, v = 1,...,M-N G . (55) 

The fact that the elementary excitations are related to the trapped pair energies can 
be readily seen in Fig. where we show the low-lying excited states for the same system 
as in Fig. 0] We see clear evidence in the figure of a phase transition, which takes place 
at roughly g ~ 0.3. For lower values of g, the states of the system are classified by the 
single-particle configurations. After the transition, this is no longer the case. We claim 
that this transition is from a normal (essentially uncorrelated) system (at small g) to a 
strongly correlated superconducting system (for large g). That crossing s take place aroun d 



the transition region is a unique characteristic of an integrable model ([Arias et all |2003). 
Observation of level repulsion in the spectrum would immediately signal non-integrability. 

In the extreme superconducting limit (g — > oo), the states with the same number N G of 
excitations are degenerate. Moreover, the slope of the excitation energies in that limit is 
given by N G - 

The degeneracies of the states in the extreme superconducting limit cIl,m,Ng have been 
obtained by Gaudin using the fact that the Richardson model maps onto the Gaudin magnet 
in this limit, and are given by 

d L ,M,N G = Cn g - C N G -n < iV G < M , (56) 
where is the combinatorial number of N permutations of L numbers. They satisfy the 
sum rule = Y,n g =o ^L,M,N a j so ^ na ^ ^ e sum °f degeneracies is the total degeneracy. 

In general, the practical way to solve the Richardson equations is to start with a given 
configuration at g — and to let the system evolve with increasing g. Hence the problem 
is to find for each initial state the number of roots No that remain finite in the g — ► oo 
limit. This is a highly non-trivial problem as it connects the two extreme cases of g = and 
g = oo. The algorithm that relates each unperturbed configuration to N G has been worked 
out by iRoman et al\ (J2003J) in terms of Young diagrams. We will not give further details 
here, but show in Fig. El a particular example with N G = 3 to illustrate the complexity of 
the evolution of the real part of the pair energies as a function of g. For sufficiently large 
values of g and after a complicated pattern of fusion and splitting of roots (2 real roots <-> 
1 complex root), the final result of N G = 3 emerges. 
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These results show the non-trivial nature of the elementary excitations of the pairing 
model, as exemplified by their non-trivial counting. As noted above, the elementary excita- 
tions satisfy an effective Gaudin equation. They also satisfy a dispersion relation similar to 
that of Bogolioubov quasipar ticles. For the s e reas ons, this new type of elementary excitation 



has been called qaudin os by 



partially addressed by 



Roman et al. 



Yuzbashvan et al. 



(2003). An interesting problem, which has been 



(2003), would be to analyze in detail the relation 



between gaudinos and Bogolioubov quasiparticles for large systems. 



VI. APPLICATIONS 



A. Ultrasmall superconducting grains 



Anderson 



(1959) made the conjecture that superconductivity must disappear for metallic 
grains when the mean level spacing d, which is inversely proportional to the volume, is 
of the order of the superconducting (SC) gap in bulk, A. A simple argument supporting 
this conjecture is that the ratio A/d measures the number of electronic levels involved in 
the formation of Cooper pairs, so that when A/d < 1 there are no active levels accessible 
to build pair correlations. Apart from some theoretical studies, this conjecture remained 
la rgely unexplored until the re c ent fa brication of ultrasmall metallic grains. 



Ralph. Black and Tinkhaml ()1997f ) (RBT), in a series of experiments, studied the super- 



conducting properties of ultrasmall Aluminium grains at the nanoscale. These grains have 
radii ~ 4-5 nm, mean level spacings d ~ 0.45 mev, Debye energies up ~ 34 mev and charging 
energies Ec ~ 46 mev. Since the bulk gap of Al is A ~ 0.38 mev, this satisfies Anderson's 
condition, d > A, for the possible disappearance of superconductivity. Moreover the large 
charging energy Ec implies that these grains have a fixed number of electrons, while the 
Debye frequency gives an estimate of the number of energy levels involved in pairing, namely 
Q = 2u£>/d ~ 150, which is rather small. Among other things, RBT found an interesting 
parity effect, similar to what is known to occur in atomic nuclei, whereby grains with an 
even number of electrons display properties associated with a SC gap, while the odd grains 
show gapless behavior. 

These experimental findings produced a burst of theoretical activity focused on the study 
of the pairing Hamiltonian [Eq. (jjjj) ] with equally spaced levels, i.e., Ej = jd. Many 
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different approa ches were used to stud y this model including: i) BCS approximation pro- 



j ected on parity ( 



von Delft et al 



Braun and von Pelf 



Mastellone et 



al 



J 19981) . hi 



19961 ) . ii) number- conserving BCS approximation (PBCS) 



Lanzcos diagonalization with up to Q = 23 energy levels 
1998). iv) Perturbative Renormalization Group combined with small di- 



agonalization ((Berger and Halperin . 



evels 



1998 ), v) the Density Mat rix Renormalization Group 



Dukelskv and Sierra 



1999), etc. [For a review on this 



(DMRG) with up to Q = 400 

topic, see Ivon Delft and Ralph! (|200l[ ).] Following this flurry of work, it was suddenly real- 
ized that the pairing model under investigation had in fact been solved exactly a la Bethe by 
Richardson long before. This came as a surprise and led to a posteriori confirmation of the 
results obtained by the "exact" numerical methods, namely Lanzcos Diagonalization and 
the DMRG. Moreover, the rediscovery of the Richardson solution produced other important 
new developments, including generalization of its solution, new insight from the point of 
view of integrable vertex mod els, connectio n with Conformal Field Theory, Chern-Simons 
Theory, etc. [For a review, see 



Sierra 



mm]. 

Returning to the application of the Richardson solution to ultrasmall superconducting 
grains, we shall focus on two quantitie s, the condensation energy and the Matveev-Larkin 



parameter ((Matveev and Larkin . 



19971 ). The condensation energy is the difference between 



the ground state energy of the pairing Hamiltonian and the energy of the Fermi state (FS), 
namely the Slater determinant obtained by simply filling all levels up to the Fermi surface. 
It is given by 



E° = Ef s - (FS\H BCS \FS) , (57) 

where 6 = for even-parity grains and b = 1 for those with odd parity. In the BCS solution, 
appropriate when the number of electrons N is very large, the leading-order behavior of 
is given by —A 2 / (2d), where A is the BCS gap in bulk and d scales as 1/N, suggesting 
that E^ is independent of the parity h of the system. However, an odd ultrasmall grain 
has a single electron occupying the level nearest to the Fermi energy. One can easily show 
that this electron decouples from the dynamics of the pairing Hamiltonian, since the pairing 
interaction only scatters pairs from energy levels that are doubly occupied to those that are 
empty. Hence the single electron only contributes through its free energy. Furthermore, 
since there is one less active level at the Fermi energy, it is harder for the pairing interaction 
to overcome the gap and the total energy thus increases. This is the physical origin of the 
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parity effect in superconducting grains. 

Recall that the BCS gap in bulk is given by A = dN / sinh(l/g), with g = 0.224 for 
Aluminum grains. The BCS result is obtained by solving the gap equation with a finite 
number of energy levels N. For even grains, there is a critical value of the ratio d^/A = 3.53, 
above which there is no solution to the gap equation. If the grains are odd, the singly 
occupied ( "blocked" ) level must be eliminated from the Hamiltonian, and the critical ratio 
becomes d\j A = 0.89. The fact that this is smaller than the even critical ratio indicates that 
odd grains are less superconducting than even grains. Thus BCS provides an explanation 
of the parity effect observed by RBT. At the same time, it suggests the existence of an 
abrupt crossover between the superconducting regime and the normal state, as conjectured 
originally by Anderson. 

However, the BCS ansatz does not have a definite number of particles, which it only fixes 
on average. Though irrelevant for a macroscopic sample, this can be important for systems 
with a small number of particles where fluctuations in the phase of the superconducting 



Braun and von D elft 



order parameter may destroy the superconductivity. For this reason 
( 19981 ) considered the PBCS state, which includes number projection and thus does not suffer 
from this limitation. There are several important differences between the results obtained 
with number projection (PBCS) and without (BCS). Firstly, the condensation energies 
from the PBCS ansatz are much lower than those from a corresponding BCS treatment. 
Secondly, the sharp transition between the SC regime and the fluctuation-dominated (FD) 
regime that arises in BCS is smoothed out by PBCS. Nevertheless, some BCS features 
survive the inclusion of number projection, particularly for odd grains. Lastly, in contrast 
to BCS, there is always a solution to the PBCS equations. 

In the upper panel of Fig. [7[ we compare the results of the condensation energy for even 
and odd grains as a function of the grain size calculated in the PBCS approximation and 
exactly. The exact solution shows a completely smooth SC/FD transition, although one can 
still talk about two asymptotic regimes which match near the level spacing for which the 
Anderson criterion d/A ~ 1 is satisfied. 

Another characterization of the parity effect is in terms of the gap parameter, which 
measures the difference between the GS energy of an odd grain and the mean energy of the 
neighboring even grains obtained by adding and removing one electron, 
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Aml = E 1 (N) - X - (E (N + 1) + E (N - 1)) 



(58) 



The lower panel in Fig. [7| compares the value of A ML / A computed at the level of PBCS 
approximation and exactly. Both curves show a minimum in this quantity as a function 
of d/A. This latter feature was first conjectured by Matveev and Larkin (1997) which 
is why the a ssociated gap is called t he Matveev-Larkin parameter. It was subsequently 



confir med by 



by 



Mastellone e.t al 



(1998) using the Lanczos method, by iBerger and Halperin 



(1998) u sing the Perturbative Renormalization Group combined with small diagonali z ation . 



Braun and von Delft 



( 19981 ) using the PBCS method, and by 



Dukelskv and Sierra 



( 20001 ) 



using the DMRG method. The shape of the exact curve, which is identical to that obtained 
using the DMRG, is rather smooth as compared with the PBCS method. This can be 
interpreted as a suppression of the even-odd parity effect. 

Richardson's exact solution of the PM can also be used to study the interplay of random- 
ness and interactions in a non-trivial model, by examining the effect of level statistics on 
the SC/FD crossover, as reflected for example in the location of the critical level sp acing. 



There was an earlier study of the latter question by iSmith and Ambegaokarl (jl99fi[ ) using 
the BCS approach, who concluded that randomness enhances pairing correlations. More 
specifically, they compared the results for a pairing model with a random spacing of levels 
(distributed according to a gaussian orthogonal ensemble) with one having uniform spacings. 
What they showed is that for both models the BCS theory gives rise to an abrupt SC/FD 
crossover, that the random-spacing Hamiltonian produces a lower correlation energy E% 
than the corresponding uniform-spacing model, and that the average value of the critical 
level spacing in the model based on random splittings is larger than the corresponding value 
for the uniform-spacing model. As noted earlier, however, the mean-field BCS theory pro- 
duces an abrupt vanishing of E^ that is not present in more sophisticated treatments. This 
raises the question of whether the conclusions they found regarding the role of randomness 
may be an artifact of the BCS approach. 

Indeed, the exact results shown in Fig. IHlfor random levels show that the SC/FD crossover 
is as smooth as for the case of uniformly-spaced levels. This means, remarkably, that even 
in the presence of randomness, pairing correlations never vanish, no matter how large d/A 
becomes. Quite the contrary, the randomness- induced lowering of E° is found to be strongest 
in the FD regime. 
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B. A pictorial representation of pairing in a two-dimensional lattice 



We next apply the electrostatic ana 



tice with a residual pairing interaction f D ukelskv et al 



ogy to a mode 



of electrons in a two-dimensional lat- 



2001). Assuming a nearest-neighbors 
hopping term, the single-electron energies in momentum space are Ek = —2 (cos k x + cos k y ), 
with k a = 2irn cr /P and —P/2 < n a < P/2. In this expression, a = x,y and P is the 
number of sites on each side of the square lattice. In the numerical example that follows, 
we consider a 6 x 6 square lattice at half filling (M = 18) with a constant and attractive 
pairing Hamiltonian for which £k = Vk- Table III shows the corresponding information on 
the positions and charges of the orbitons in the subspace of seniority-zero states. 

Table III: Positions and charges of the orbitons for an attractive 2D pairing model. 
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Figure |H1 shows the equilibrium positions of the pairons associated with the ground-state 
solution for three values of g. The orbitons are represented by open circles with radii 
proportional to their charges, and the pairons are represented by solid circles. For each 
pairon in the figure, we draw a line connecting it to the one that is closest to it. 

In the limit of weak pairing {g = —0.040), 5 pairons are very near the half-filled orbiton 
that is located at the Fermi energy (E = 0) and has charge —5. The other 13 pairons are 
distributed close to the lowest four orbitons, consistent with the Pauli principle. There is one 
near the lowest and four near each of the next three. This part of the figure suggests a picture 
whereby for weak pairing the pairons organize themselves as artificial atoms around their 
corresponding orbitons. As g increases, the pairons repel, causing the atoms to expand. For 
g = —0.064 the two orbitons closest to the Fermi energy have lost their pairons which have 
now linked up with the five that are near to the third orbiton. [We can make this linkage more 
precise by drawing lines that connect each pairon with its nearest neighbor.] We refer to this 
grouping of 13 pairons as a cluster, as it arises from the pairons that originally comprised 
three atoms. The remaining 5 pairons are still attached to their orbitons, as artificial atoms. 
By g = —0.130, the cluster has grown to the point that all 18 pairons are trapped in it. 
We claim that this delocalization effect in the classical problem, from independent atoms to a 
collective cluster, is a pictorial reflection of the transition from a normal to a superconducting 
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system in the quantum problem. In the extreme superconducting limit, as reflected in the 
figure by the g = —0.130 panel, all pairons are behaving collectively and have lost their 
memory of the orbitons from which they arose. In the corresponding quantum problem, 
there are a set of Cooper pairs that likewise behave collectively and have lost their connection 
to specific single-particle orbits. 



Similar results have been reported by 



Dukelskv et al 



(2002) for the problem of pairing 



between alike nucleons in atomic nuclei. The analysis was for two isotopes of Sn, including 
the isotope U4 Sn briefly alluded to in Sect. III. There too the superconducting phase tran- 
sition was seen to be associated with a transition from isolated atoms to a cluster in the 
analogous electrostatic picture. Furthermore, there too the transition to full superconduc- 
tivity was seen to develop in steps, depending on the single-particle levels that play a role 
in producing the pair correlations and their ene rgy hierarchy. 

The analysis of pairing in nuclei reported by iDukelskv et a 
interaction. Of course, this is just an approximation to the true nuclear interaction in the 
J 71 " = + channel. Nevertheless, we expect that the general features of the transition to su- 
perconductivity should be the same even for a more realistic pairing interaction. It is in fact 
possible to build greater flexibility into the nuclear structure analysis, while still preserving 
the electrostatic analogy, by considering more general exactly-solvable Hamiltonians of the 
rational family. 



( 20021 ) assumed a pure PM 



C. Electrostatic image of interacting boson models 

As an example of the use of the electrostatic mapping for a finite b oson system, we now 



discu ss the phenomenological Interacting Boson Model (IBM) of nuclei (jlachello and Arima . 



198CM ) . The IBM captures the collective dynamics of nuclear systems by representing corre- 
lated pairs of nucleons with angular momentum L by ideal bosons with the same angular 
momentum. In its simplest version, known as IBM1, there is no distinction between protons 
and neutrons and only angular momentum L = (s) and L = 2 (d) bosons are retained. 
We will use the electrostatic image to study the properties of a second-order quantum 
phase transition that arises in the IBM1 from a vibrational system with U(5) symmetry 
to a gamma-unstable deformed system with 0(6) symmetry. This phase transition can be 
modelled by the one-parameter IBM1 Hamiltonian 
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H= n d + |ptp , (59) 

where n d = c^c^, pt = s t s t _ ^ (— 1) M <$d_^ creates a boson with angular mo- 
mentum L = 0, dl creates a boson with angular momentum L = 2 and z-projection /i 
(—2 < // < 2), and x is the ratio of the pairing strength g to the single-particle splitting 
Ed — e s . The parameter x can be varied from x = (the £7(5) limit) to x = oo (the 0(6) 
limit). Equation (J59j) is an example of an exactly-solvable repulsive pairing Hamiltonian, 
and the second-order natu re of the phase tra nsition it describes has been recently attributed 



to quantum integrability ([Arias et all 120031 ) . 

The electrostatic problem that corresponds to this quantum boson model consists of two 
orbitons with positive charges q s = 1/4 and q d = 5/4 (see Sect. III). Both the s and d 
orbitons are located on the real axis, with the s located at position 0.0 and the d at position 
2.0. There are M pairons with positive unit charge that interact with the orbitons and 
with one another and that feel an external electric field pointing downwards with strength 
1/x. In the ground-state configuration, the pairons are constrained to move between the 
two orbitons. 

Figure El shows the pairon positions for a system of 10 bosons as a function of the 
control parameter x. For x close to 0, corresponding to weak repulsive pairing, there is 
a very strong electric field, which compresses the pairons very close to the s orbiton. As 
x increases, the electric field decreases and the Coulomb repulsion among all the charges 
begins to counteract the effects of the external field. As a consequence, the pairons gradually 
expand along the energy interval between the s and the d orbitons. The phase transition 
between the vibrational U(5) system and the gamma-unstable 0(6) system, clearly depicted 
in the classical electrostatic problem, arises when the Coulomb repulsion and the external 
electric field balance one another. Prior to the phase transition, the quantum system is 
primarily an s boson condensate, with perturbative contributions from d bosons. Following 
the phase transition, it is a fragmented condensate m ixing the s and d bosons. 



A similar analysis by 



Dukelskv and Pittell ()200lh has also been carried out for a system 



with many even angular momentum boson degrees of freedom, not just the s and d. The 
purpose of that analysis was to better understand why the IBM1, with just s and d bosons, 
works so well in describing the collective properties of nuclei. Recognizing that bosons in this 
model are an ideal representation of the lowest fermion pairs of identical nucleons and that 
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there are not just + and 2 + pairs, a natural question to ask is: Why can we ignore the higher 
angular momentum pairs/bosons when dealing with nuclear collective properties? Part of 
the answer is contained in the dominant quadrupole-quadrupole neutro n-proton interaction, 



Dukelskv and Pittel 



which is known to favor the lowest + and 2 + pair degrees of freedom 
(2001) suggested another mechanism, based on an analysis of a generalized boson model con- 
taining all even angular momenta up to some maximum and interacting via a repulsive boson 
pairing interaction. The latter is a means of simulating the repulsive interaction between 
bosons that arises due to the Pauli principle between the fermion (nucleon) constituents of 
which they are comprised. Using the Richardson solution of this model, they showed that 
a repulsive boson pairing interaction can only correlate two boson degrees of freedom, and 
that these should be the lowest two, the s and th e d. More recently, this re sult has been 
interpreted by means of the electrostatic mapping ((Pittel and Dukelskvl . 120031 ). Even in the 
presence of many boson degrees of freedom and thus many orbitons, the collective pairons 
are always confined to lie between the lowest two, i.e., between the s and the d. 



D. Application to a boson system confined by a harmonic oscillator trap 

We now consider the problem of a set of bosons confined to a harmonic oscillator trap and 
subject to a boson pairing interaction. We claim that such a Hamiltonian cannot realistically 
describe the physics of a confined boson system, for the following reason. Looking back at 
the commutators of the pair operators A\ given in Eq. (J2J), we see that they are normalized to 
the square root of the degeneracy fli of the level I. Thus, the PM Hamiltonian of Eq. (JIJ has 
a pairing matrix element proportional to y/QiQi>. In a three-dimensional harmonic confining 
potential, these degeneracies are in turn proportional to I 2 , where / plays the role of the 
principal quantum number and the summation in the pair operators of Eq. (JIJ) now include 
both the orbital and the magnetic quantum numbers. On the other hand, the single-boson 
energies E\ for such a confining potential are linear in I. Thus, a boson pairing interaction in 
the presence of an oscillator confining trap would have the net effect of scattering boson pairs 
to high-lying levels with greater probability than to low-lying levels, producing unphysical 
occupation numbers. 

To numerically test this conjecture, we have solved the Richardson equations [Eq. (|34j) ] 
for a system of 1000 bosons (M = 500) trapped in a three-dimensional harmonic oscillator 
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(ftj = (l + l)(l + 2) 12 and e, = hul l + 3/2))with a cutoff at 101/2^ (L = 50 single boson 
levels). Following iRichardsonl ([19681 ). the occupation numbers can be calculated as 



From Eqs. (@J) and (|34|) . a set of M coupled nonlinear equations in terms of M new 
unknowns are obta ined, which when solved give the L occupation numbers. For details of 



the derivation , see 



Richardson! (jl968h 



In Fig. we show the occupation numbers versus the single-boson energies in units of 
hut for a pairing strength of g = —0.0025. We have excluded the occupation of the I = 
condensed boson state from the figure, since it lies outside the scale of the figure. The overall 
depletion is 0.21, which gives an occupation of the I = state of Uq = 790 bosons. The figure 
clearly shows that the depletion is unphysically dominated by the high-lying (i.e., high-/) 
levels, due to the nature of the pair coupling matrix elements discussed above. 

We can use the freedom we have in choosing the parameters rji entering in the definition 
of the R operators to obtain a more physical exactly-solvable model. In order to cancel the 
unphysical dependence of the pair coupling matrix elements on the degeneracies, we choose 
the 77's so that rji = (si) 3 . Then, the new Hamiltonian, which is given by the associated 
linear combination of R's, will be 

H = 2J2 eiRi = C + J2 mi + E Vw [A\A V - mm] , (61) 

l l lj4> 

where 

1 1 4 w 

si = si- E VwSlv, 



Vu> = £ -^—TV^- ■ (62) 



9 1_ 

2 ef + ef, + EiE V 

Taking into account that E[ is proportional to I, the two-body matrix elements in Eq. (|62[) 
cancel the dependence on the degeneracies in the effective pair coupling matrix elements. 
Thus the above Hamiltonian should be more appropriate when modelling a harmonically- 
confined boson system with a pairing-like interaction. 
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The interaction in Eq. (|62|) has the nice feature that its two-body matrix elements 
decrease with the number of shells, as one would expect in general. It has the particular 
property that the interactions of the pair- and density-fluctuations are strictly the same 
but opposite in sign. The energy eigenvalues of this Hamiltonian can be obtained from the 
eigenvalues r ; of the associated Ri operators as E = 2 J2i £ i r u with the end result being 

E = W * A - 7 E VwWv -2gJ2 -f^- . (63) 

Z l ^ tyl' lp A£ i 

Note that the first two terms in the energy eigenvalues of Eq. (p)3"j) exactly cancel the constant 
term C in Eq. (jfiT]). 

We have solved the Richardson equations for the Hamiltonian of Eq. ()62jl . i.e., with 
Tji = (^/) 3 , for the same system considered above, namely M =500 and L =50. In Fig. El 
we show the occupation numbers for g = —1.0, for which the overall depletion factor is 0.245. 
They display a reasonable physical pattern, with the occupancies decreasing monotonically 
with increasing single-boson energy. A comparison between the exact results and such 
approximations as Hartree-Fock-Bogolyubov theory or its number-conserving variants will 
be the subject of future work. 



When the same modified form for the Hamilto nian, but with repulsiv e 



pairing, was con- 



sidered, a highly unexpected feature was found ((Dukelskv and Schuckl . l200lh . Figure [T3] 
shows the occupation numbers of the first and second levels versus the scaled pairing in- 
teraction x = 2Mg/huj for the same system as above (M = 500, L = 50). At the critical 
value of the scaled interaction, x c — 1, the normal ground-state boson condensate suddenly 
changes into a new phase in which the bosons occupy the I = and the / = 1 levels, with 
the occupation of all other levels negligible. 

This new phase is characterized in the large- N limit by having two macroscopically oc- 
cupied stat es, thus representing a fragmente d condensate. It is commonly accepted since 



the work of iNozieres and Saint Jamed (|l982t ) that for confined systems fragmentation can- 
not occur in systems of scalar bosons with repulsive interactions. This might be the first 
example of fragmentation in a confined boson system. 
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VII. SUMMARY AND OUTLOOK 



In this Colloquium, we have reviewed recent efforts to develop exactly-solvable models of 
the Richardson-Gaudin type and have discussed how these models have been used to provide 
valuable insight into the physics of systems with strong pair correlations. We began with a 
brief review of Richardson's original treatment of the pairing model and Gaudin's related 
treatment of the Gaudin magnet and then showed how they could be generalized to several 
classes of exactly-solvable models that still preserve a pairing-like structure. 

A very attractive feature of these models is that one can establish an exact analogy 
between the associated quantum many-body problem and the classical physics of a two- 
dimensional electrostatic system. This feature, which was originally appreciated by both 
Richardson and Gaudin, has now been generalized to all such exactly-solvable models. 

The solution to a classical problem is typically amenable to simple geometrical interpre- 
tation, which can then be used to give an alternative perspective on the quantum model 
from which it derives. This has been used in several examples we have discussed, both for 
boson and fermion systems. It provides a new perspective on how superconductivity arises 
in fermion models and also an interesting new perspective on a second-order quantum phase 
transition that arises in the nuclear Interacting Boson Model. 

Another important outcome of the classical electrostatic analogy is that it facilitates 
a treatment of these models in the thermodynamic limit. Not surprisingly, it shows that 
BCS theory is indeed the correct large- A limit of the fermion pairing model with attractive 
interactions. 

The generalization of Richardson and Gaudin's models was first discussed in the context 
of the physics of small metallic grains, as a means of obtaining exact solutions to the BCS 
Hamiltonian appropriate to such systems in cases where numerical diagonalization was out 
of the question. As discussed in some detail in this Colloquium, the existence of a method 
of exact solution for this model provides important insight into the detailed physics present 
when the size of the system gets small enough so that superconductivity disappears. Neither 
BCS theory nor number-projected BCS theory can capture the physics of the phase transition 
acceptably. 

The new families of exactly-solvable models that we have discussed are based on the 
algebras SU(2) (for fermion systems) or SU(1,1) (for boson systems). These two algebras 
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have a pseudo-spin representation in terms of fermion pairs or boson pairs, respectively. All 
applications to date and all that we have therefore discussed are based on these represen- 
tations of the associated algebras. At the same time, the SU(2) group has other possible 
representations in terms of spins operators or two-level atoms which have not yet been 
exploited. 

Perhaps, the most important feature of the exactly-solvable R-G models is the enormous 
freedom within an integrable family. For a given set of L single-particle levels (or orbits), 
we can define an exactly-solvable Hamiltonian in terms of L + 1 r/j and g parameters that 
enter in the definition of the associated quantum invariants [Eqs. 129|H0|) ]. and an additional 
set of L parameters £j that define a Hamiltonian as a linear combination of the quantum 
invariants. The models therefore contain 2L + 1 free parameters, which allows enormous 
flexibility i n constructing a genera l pairing-like interaction tailored to the physical problem 



of interest fD ukelskv et al 



2003). In contrast, most other exactly-solvable models have 



either no free parameters (the Heisenberg model), one free parameter (the XXZ model, the 
Hubbard model, or the Elliott model), or just a few free parameters (the three dynamical 
symmetry limits of the nuclear Interacting Boson Model). We reported here one particu- 
larly interesting use of this flexibility, namely to model the physics of bosons confined to a 
harmonic trap. A pure pairing interaction has the anomalous feature that it scatters pairs 
of bosons preferentially to high-energy states. By exploiting the flexibility of the general- 
ized R-G models, it was possible to find a more physically-meaningful Hamiltonian for such 
systems, which nevertheless was still exactly solvable. Analysis of this new model led to a 
suggestion that confined boson systems can exist in a fragmented state, contrary to prior 
expectations. 

All the application reported discussed in this Colloquium made use of generalized R-G 
models of the so-called rational class. Other interesting applications within this class of 
models should certainly be sought. 

A potentially interesting example concerns the properties of the pairing phase transition in 
finite Fermi systems at finite te mperature. There has already be en important work reported 
on this topic (see, for example, iDean and Hiorth- Jensen! (J2003J) and reference therein). As 
noted earlier, the exact solution of the pairing model can be derived from the Algebraic 
Bethe ansatz. It is natural, therefore, to study the Thermodynamic Bethe ansatz, which 
provides the exact finite temperature description of short-range ID integrable models, to 
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see whether it could be extended for application to integrable pairing models. If so, this 
would allow for an exact treatment of the finite temperature properties of such systems as 
nuclei, ultrasmall superconducting grains and degenerate Fermi g ases. 



Another important topic that has recently received attention ((Dean and Hi orth- Jensen . 



20031 1 is the study of the low-energy properties of quantum many-body systems with random 



interactions and in particular with random pairing interactions. Though random pairing 
interactions in general have chaotic properties, there is an important subset of integrable 
pairing hamiltonians that have a large number of free parameters which can be chosen 
randomly. Some of the physical consequences of randomly chosen pairing interactions in 



the general chaotic regime and within the integrable subset o 



reported recently by 



Volva et al 



( 2002h and 



Relano et al. 



pairing models have been 



(2004), respectively. 



We also expect interesting applications to ensue for the t rigonometric a nd hyperbolic 
models. As one example, the trigonometric model was used by Gaudin (jl976 ) to derive the 
limit of an exactly-solvable model of the interaction of a two-level atomic system with an 
external oscillator field. This kind of model could lead to a generalization of the celebrated 
Jaynes-Cummings model. Using techniques similar to the Algebraic Bethe ansatz derivation 
of the Richardson model, exactly-solvable models for boson at omic-molecule syste ms and 
two coupled Bose-Einstein condensates have recently been found ([ Links et all 120031 1. Based 
on these arguments, we believe that R-G models could have a promising future in Quantum 
Optics and in the study of dilute Fermi and Bose gases. 

Another area of great interest is the generalization of R-G models based on the SU(2) or 
SU(T,1) algebras to larger algebras. Some work has already been reported along these lines 



f Asorev et al 



20021 ) . A complete solution for models with 0(5) or SU(4) symmetries could 



lead to interesting applications for nuclear systems with N m Z, where it is important to 



explicitly include the isospin degree of 



prope rties of high T r supercondu ctors (jWu et al 



works 



Ric hardson! (|1366b, 



reedom. 



t could also provide useful insight into the 



200 



3; 



Zhang I . Il997f l . Though in previous 



1967) proposed an e xact solution for p airing Hamiltonians with 



these two group symmetries, recent work by IPan and Draaverl (1200211 indicat e s tha t his 



solution is invalid for systems with more than two pairs. Moreover, 



Guan et al. 



Links et all (1200211 and 



( 2OO2I ) have found different solutions for the same problems. More work is clearly 
required along these lines. 

In closing, the use of exactly-solvable Richardson-Gaudin models has already provided 
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significant new and important insight into the properties of many diverse quantum systems, 
ranging from atomic nuclei to electronic systems in condensed matter. We are optimistic 
that many more exciting applications are still to come. 
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FIG. 1 Two-dimensional representation of the pairon positions in 114 Sn for a pure PM Hamiltonian 
with single-particle energies as given in Table II and with g = —0.02 MeV. 
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FIG. 2 Evolution of the real and imaginary parts of E^{g), in units of the mean level spacing 
d = uj/L, for the equally-spaced model with M = L/2 = 8, as a function of the coupling constant 
g. For convenience, the energy levels are chosen in this figure as €j = 2j. The dashed curves in 
the lower half are an average of the two neighboring real energies. Following the criti cal point, this 
turns into the real part of the energy of the resulting complex conjugate pair. From 



Roman et al. 



(2002). 
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FIG. 3 Plot of the roots for the equally-spaced pairing model in the complex £-plane. The 
discrete symbols de note the numerical values for M = 100. All energies are in units of u. From 



Roman et al 



(2002) 
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FIG. 4 Real part of in units of the mean level spacing d for the equally-space d pairing model 
with M = L/2 = 20 pairs and Nq = 0, 1, 2 excitations. From 



Roman et al. 



(2003). 
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FIG. 5 Lowest 44 excitation ener gies E P . xr = E — Eg g i n units of the mean level spacing d for 



M = 20 pairs at half filling. From lRoman et al. 



(2003). 
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FIG. 6 An illustration of the algorithm for determining the number of roots Nq, for a case with 
Nq = 3. a) Top: Initial state at g = 0, where each solid circle (•) denotes a level occupied by a 
pair, each empty cicle (o) denotes an empty level, and the vertical dashed line with an F near it 
denotes the position of the Ferm i surface. Bottom: The associated Young diagram, b) Real part 



of E„. From 



Roman et al. 



(2003). 
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FIG. 7 Calculated results for ultrasmall superconducting grains as a function of the grain size. 
The upper panel gives the condensation energies £f for even (b = 0) and odd (b = 1) grains 



calculated using PBCS and exact wav e functions. The lower 



A^, obtained in PBCS calculations ( Br aun and von Delft . 



pane l gives the Matveev-Larkin gap, 



1998) and in exact calculations. The 



quantity A that scales the condensation energy, the Matveev-Larkin gap and the mean level spacing 
d is the BCS gap in bulk. 
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FIG. 8 Exact even and odd condensation energies Eg for uniform equally-spaced levels (dashed 
line), and the ensemble average energies (Eg) for random-spaced levels (solid line). The height of 
the fluctuation bars gives the variances 5 Eg. The results are plotted as a function of the mean 



level spacing, d, scaled by the BCS gap in bulk, A. From 



Sierra et al. 



(2000). 
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FIG. 9 Two-dimensional representation of the positions of the orbitons and pairons corresponding 
to a 6 x 6 lattice at half filling, for three selected values of g. The orbitons are represented by open 
circles and the pairons by solid circles, with radii proportional to their charges. 
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FIG. 10 Evolution of the pairon positions as a function of the scaled strength parameter x for 
a model with 10 s and d bosons subject to a Hamiltonian with linear single-boson energies and 
a repulsive boson pairing interaction. The circles at x = denote the positions of the s and d 
orbitons. 




FIG. 11 Occupation numbers for 1000 bosons confined in 50 harmonic oscillator shells and inter- 
acting via a pure pairing interaction with strength g = —.0025. The occupation of the 1=0 state is 
off scale and thus not shown. 
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FIG. 12 Occupation numbers for 1000 bosons confined in 50 harmonic oscillator shells and inter- 
acting via a renormalized pairing interaction [see Eq. (|6'2|)] with g = —1.0. 
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FIG. 13 Occupation numbers no and m for 1000 bosons confined to 50 harmonic oscillator shells 
and interacting via a repulsive renormalized pairing interaction as a function of the scaled strength 
parameter x. 
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